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Abstract 

We present an efficient method for building equilibrium multi-component galaxies 
with non-spherical haloes and bulges. The gist of our approach is to tailor the 
velocity ellipsoid directly to the geometry of the mass distribution. Thus we avoid 
computing the anisotropic velocity dispersions which leads to large savings in the 
cpu budget. The computational time of the algorithm for triaxial equilibria scales 
linearly with the number of particles, A^. The approximate solution to the velocity 
field causes structural relaxation: tests with = 50, 000 show that fluctuations 
of the inertia tensor (not exceeding the 10% level) disappear within one half of a 
revolution at the half-mass radius. At later times equilibrium properties settle to 
values close to those sought from the initial conditions. 

A disc component is then added as described by Hernquist (1993). Incorporating 
the above algorithm to his code buildgal , test runs show that the stability of the 
disc against vertical heating is not substantially modified by using our method. 

The code, MaGalie, is made generally available^. 

Key words: Galaxies, dynamics; halo; dark matter; numerical: method 
PACS: 95.35+d, 95.75.Pq, 98.52-b 



Email : cmb@ari.uni-heidelberg.de 
^ Email : pavel@astrophysik.uni-kiel.de 
^ Email : jorpega@ari.uni-heidelberg.de 



New Astronomy, in press 



1 February 2008 



1 INTRODUCTION 



The visible Milky Way and other disk galaxies may be contained in extended 
dark-matter haloes of largely unknown morphology. Spherical haloes are typ- 
ically assumed, but there are indications that this matter component may be 
significantly fiattened, or even triaxial (see Zaritsky 1998, Sackett 1999, Rusin 
& Tegmark 2000). An important indicator is the non-isotropic distribution of 
satellites around disc galaxies, or the Holmberg effect (see Holmberg 1969). 
This may be a result of dynamical friction of non-spherically symmetric dark 
haloes acting on the satellites (Penarrubia, Boily & Kroupa 2000). 



Equilibria of self-gravitating galaxies are achieved for a given distribution func- 
tion (df) based on three integrals of the motion. This approach faces a dif- 
ficulty in that only the first two classical integrals (total energy E and one 
component of angular momentum J) are known, the final integral taking an 
assumed form (see e.g. Lupton, Gunn & Griffin 1987; Einsel & Spurzem 1999). 
A greater difficulty still is that galaxies are heterogeneous systems of several 
components, and hence are hardly represented adequately with a unique df. 
A self-consistent equilibrium may be constructed by resorting to numerical 
integration of the equations of motion. Barnes (1988) introduced an iterative 
scheme by which particles relax in a fixed potential representing a galaxy 
component. This component is then given a particle representation in near- 
equilibrium. 



Hernquist (1990, 1993) has introduced a by-now standard scheme which solves 
moments of the coUisionless Boltzmann equation using tailored distribution 
functions. Three- component galaxies so constructed have the advantage over 
the Barnes approach of closer agreement with the sought equilibrium. In prac- 
tice his algorithm is highly efficient when embedding a galactic disc in spheri- 
cal components. (Benchmarking examples will follow in subsequent sections.) 
However we have found that for the purpose of constructing axisymmetric 
(non-spherical) bulges and haloes, the scheme requires significantly more cpu- 
time, mainly because one must solve numerically for the components of the 
velocity dispersion locally in three directions. This introduces an A^^ scaling 
of the computational time. By contrast, in spherical symmetry the velocity 
dispersion is known immediately from a one-dimensional integral (cf. [17] be- 
low). The algorithm of Hernquist (1993) remains efficient for up to tens of 
thousand of particles, however the computational time becomes prohibitive 
for larger numbers. Given that particle numbers adequate to galactic systems 
by far out-stretch available computer power, algorithmic problems such as this 
require solutions. 



2 



In this contribution we by-pass the Boltzmann equations approach, by noting 
that galactic components follow a well-defined hierarchy of sizes and masses. 
This suggests a perturbative treatment of particle potentials from the largest 
scales, down. Equilibrium solutions to Boltzmann equations in spherical sym- 
metry are first transformed to the desired morphology directly, by mapping the 
velocity ellipsoid to the density isocontours. Wc construct equilibrium multi- 
component axisymmetric galaxy models by perturbing the velocity field of 
individual components to account for the added background potential (or, em- 
bedding) . The computational time required for the whole procedure scales with 
particle number as for spherically symmetric models. Thus multi-component 
axisymmetric equilibrium configurations with millions of halo particles are 
easily constructed. In the following Section we review the problem and mo- 
tivate our approach. This is then implemented and tested using a grid code. 
Benchmarking and comparisons with Hernquist's standard buildgal follow. 
We conclude with suggestions for applications and further improvements. 



2 METHOD 

2.1 Scaling of Computational Time 

Approximate solutions based on moments of the Jeans equations have proved 
useful when constructing equilibrium galaxies. Hernquist (1993) has intro- 
duced a scheme whereby equilibrium galaxy components are first constructed 
in isolation, then zipped together by adding the gravity of the other compo- 
nents in turn. The modified binding energy of individual particles is matched 
by a similar increase in kinetic energy to preserve the structure of each com- 
ponent. 

The case when all components are spherically symmetric allows a simplifica- 
tion. Then the added binding energy is measured from the increase in mass 
inside a particle radius r. To add together two or more spherical components 
only requires to sort the particles radially. The computational time required 
by a scheme such as QUICKSORT scales no faster than the 3/2-power of particle 
number (see e.g. Press et al. 1992, §8). 

For non-spherical potentials however, the modified velocity field (and square 
velocity dispersions) requires the computation of the local potential (and its 
spatial derivatives) for each particle. The exact, direct-summation algorithm to 
do this requires a number of operations increasing in proportion to oc A^^, the 
square particle number. The particle realisation of the density itself requires 
oc operations. Thus in general we may write the total computational time 
to realise a compound equilibrium of c components as 
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Background potciitialsX 

(1) 



^ ' j=l 

y Spherical or aspherical 



where Ni is the number of particles for component i; tden, tvei the computa- 
tional times per particle to realise the density and velocity fields for respective 
components; and tcor,j the characteristic time for adjusting the velocity field of 
i*''-component particles due to non-spherical potentials. For three-component 
galaxies made of a spherical bulge (b), disc (d) and aspherical halo (h), we 
would write 



T — Nh tden,h + tvel,h + Nh {Nd + Nb) tcor,h + Nd tden,d + tvel,d + 

Nd {Nh + Nb) tcor,d + Nb tdenfi + Nb Uel,b + Nb {Nh + Nd) Uorfi ■ (2) 

Clearly massive haloes call for large particle numbers and so the leading terms 
in (2) are oc N^ and oc Nh Nd- We wish to remove these dependencies on Nh- 

In any multi-component galaxy, the largest-scale structure would be closest 
to equilibrium self-gravity. This suggests that we look for a general scheme 
whereby the relative masses and scales of components are taken into account. 
Figure 1 provides a cartoon representation of our approach. Since details of 
the small scale components are lost on the much larger ones, a correction 
to the velocity field of larger structures would only involve truncated expan- 
sions of the small-component potentials. By contrast, a light structure requires 
more accurate treatment of the surrounding mass distribution: this is shown 
as light and heavy arrows on Fig. 1. We adopted this strategy by solving for 
the self-gravitating aspherical halo through a map of the velocity field of the 
spherical solution. As we shift down the scale ladder to smaller structures, 
higher-order expansions provide accurate corrections to the velocity field of 
individual components. Thus feedbacks from disc and bulge on the halo ve- 
locity field are given a low-oder series-expansion treatment, while the velocity 
field of disc and bulge particles are adjusted with a high-order expansion of the 
halo potential. Below we take a stepwise angle to implementing this algorithm. 



2.2 Self- gravitating axisymmetric equilibria 



The basic step consists in building axisymmetric systems from a transforma- 
tion of spherically symmetric equilibria. If we consider as given the geometric 
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aspect ratio and the mass profile of tlie equilibrium sought, it is then sufficient 
to construct a velocity field to match the gravity. In cylindrical coordinates 
{R = \Jx^ + y^, 9), the orbits of stars in an axisymmetric harmonic poten- 
tial are recovered from integrating two decoupled equations of motion. In the 
case of a uniform oblate spheroid of density p, axes a,c ;a > c, eccentricity 
e = y'l — c^/a^, the potential inside the bounding volume is (Table 2-1 in 
Binney & Tremaine 1987, hereafter BT+87; Chandrasekhar 1969, §3) 



^sp{R, z) = -nGp (l{e)a^ - A^{e)R^ - A^{e)z') (3) 

with known algebraic expressions for /, Aj. The motion of a test particle is 
found from solving for 



Xi = -Vi$ = -2nGpAiXi , i = 1...3 (4) 

which admits a periodic solution in each component. Averaging the square 
velocity and position over one cycle, and substituting for Ai from Table 2.2 of 
BT+87, we obtain 



<vi> _ As <xl> _ As / ^ o "^"" ~ ^'"^ ^^^^^"^ (5) 

<Vi> Ai<Xi> Ai^ ^ sm~^{e)/e - c/a a? ' 

Prom (3) we find lines of constant potential intersecting the z — Q plane 
and the z-axis such that for these coordinates Aix\ — ^3X3. Defining the 
eccentricity of isopotential contours 





(6) 



we may relate e$ to the geometric eccentricity (e) and that of the velocity 
ellipsoid, e^, using (5) and (6) 



i 



1 - 



<vl> 
<vl> 



l-e| 



(7) 



The three quantities satisfy e$ < < e, and hence the velocity ellipsoid is 
never as flat as the mass distribution that gives rise to it. This result may be 
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derived directly from the virial theorem apphed to the system as a whole (see 
BT+87, §4.3). Therefore, for the harmonic potential, (7) applies equally to 
averaged quantities and individual orbits. 

In practice galaxies and haloes show peaked density profiles and hence the 
exact result (7) has limited application. An alternative to (7) is obtained by 
considering a star on a loop orbit in the mean potential of all the others. Since 
the matter is concentrated around the centre of gravity, a series expansion 
will be adequate when the orbit avoids the central region. Consider therefore 
the potential exterior to a homogeneous spheroid, of axes (a, c) equal to the 
mass-weighted means of the peaked distribution. This is written in a form 
analogous to (3) (BT+87, Chandrasekhar 1969), 



^,p{R, z) = -nGp- 



a'Yd 



I{a',d){a'f-Y.A,{a',d)x^^ 



In terms of (a, c) we have 



% + 7^ = 1 ; a' = sf^^TX, d = V^TX . (9) 



{a'f id) 

For oblate spheroids the solution A(a, c, r) is expressed in algebraic form: 



2 A = - a\2 - e') + ^J{a^ - c^f + 2a^e\z'' - B?) + , 

where = i?^ + z^. (A relation similar to this one also exists for prolate 
spheroids but is omitted for reasons of clarity.) The solution A > applies 
outside the bounding surface of the spheroid. For a loop orbit with r > a a 
series expansion of (8) leads to (cf. Goldstein 1980, eq. 5-88) 



$,p(it:, z) = ^ + -^{h - Ir)P2{cos9) + 0{[a/r]') (10) 

where tan 6 = z/R, Pn{x) is a Legendre polynomial, and Lj. an eigen-component 
of the inertia tensor per unit mass and we consider only the half-space 2; > in 
what follows. Integrating the equations of motion for the potential (10) leads 
to solutions of the form 
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/ ^2 . \ 

Xi{t) = Xo,iCOs[wi{t-to)] il+^-^COs'^[wi{t-to)] + ...\ (11) 



where 



Wi = Gp+ -Gp — ; W2 = -Gp — ; p = 



Thus when r 3> a > c, 'u;2 — >^ with 1/r^ since the /j's are constants. The 
orbits describe eUipses of the same eccentricity in both x and x . Note when 
h — Ir, (11) reduces to solutions for the right-hand side in (4). For such orbits 
we would therefore set 



= e . (12) 

Our solutions (7) and (12) bracket the range of possible aspect ratios the 
velocity ellipsoid may assume in equilibrium. Realistic galaxy equilibria will 
fall somewhere between strict homogeneity and point-mass distributions. We 
therefore sought a sensible interpolation for e^, in the range e$ to e to account 
for this. The number of possibilities is a priori endless. We chose for each 
particle i 



2 I / 2 2 \ 



\ 



<r2> 



g,^ 



(13) 



which was found to yield adequate equilibria. In computing e^,, we substitute 
in the above the mass- weighted average inside the spherical volume 47rr^ j/3. 
To obtain an eccentricity for the velocity ellipsoid that reflects the character 
of the gravitational field sampled by the star over one orbit, we compute the 
gravitational radius 



_ GM 

where Ei is the binding energy of the i^^ star, and M the total mass of the 
system. The velocity ellipsoid computed from (13) is the same for all stars of 
the same binding energy. We give some motivation for the interpolation (13) 
through an example. 
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To fix ideas, consider a Dehnen model for spherical galaxies. Delinen (1993) 
introduced a family of potential-density pairs which provide a suitable fit to 
bulges and cuspy ellipticals : 



M(r) = Mo 



r 




$(r) 



GMo 



X < 




1 - 



r + Tc 



r 



12-7 



) 



(7 7^2) 



(14) 



In 



r 



(7 = 2) 



Here (rc, 7) arc two parameters and the total system mass = Mo- The case 7 = 
1 corresponds to a spherical Hernquist (1990) density profile. Let Vc = 1/10, 
and we truncate the distribution at r = 1 for convenience. With Mg — 1.21.. 
the mean square radius <r'^>— 0.349..., enclosing 63% of the mass. We invoke 
the virial theorem and set Ei = — $j/2 for a given stellar orbit. Using this in 
(14) we compute Vg^i = GMo/ (— $j/2) = 2(rj + rc); the ratio of averaged values 
<rg> I <r>= 1/4. If we now squash the sphere down the z-axis to achieve an 
aspect ratio c/a — 1/2, we compute from (13) an eccentricity ~ 0.860.. , 



close to e — ^1 — / 0? — 0.866... Setting 7 = so that the central region is of 
uniform density, we find from (7) for this model = 0.832.., nearly unchanged 
from the previous value. This is because only ^ 12% of the mass now resides 
inside r = Vc- Thus unless the galaxy shows a broad uniform density region, 
the solution (13) on the mean lies close to e. 



3 Stand-alone Axisymmetric Haloes and Bulges 

We now implement the scheme described above. We begin by summarising 
the steps used in constructing spherical equilibria. Applications will centre on 
the Dehnen models (14), as they cover a range of central density gradients 
(through the index 7) from which to test the method. We concentrate on two 
aspects, namely the stability of equilibria so constructed and the range of cases 
to which this algorithm may be applied. We consider isolated self-gravitating 
haloes first; multi-component models are discussed in the next section. 

3.1 Spherically Symmetric Equilibria 

For the case of spherically symmetric mass profiles the potential at radius r is 
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$(r) = - r ^^^^ dx = - r r 47r«V(«) du (15) 

provided the run of mass M(r) increases more slowly than oc r at large dis- 
tances. The mass distribution is a given quantity of the problem. A particle 
representation of any mass profile M(r) can be obtained straightforwardly, 
first by filling a spherical volume uniformly by Monte Carlo method, then 
mapping the mass shells to the desired profile, i.e. Mh{f) = M(r), with Mh 
the uniform-density spherical profile. The solution r(M) is not known in ana- 
lytic form in general and must be tuned by interpolation when only a discrete 
representation Mj(rj) is given on input. One advantage over solving directly 
by Monte Carlo method for r(M) is that a unique realisation Mh{r) may be 
transformed to a variety of profiles. This essentially keeps any bias due to 
sfN fiuctuations the same for different M(r). This is especially important in 
systems where the effect of low-N statistics may be felt, such as in central 
galactic cusps. Interpolation errors would wipe out this advantage when the 
desired mass profile is known only on a coarse mesh {rj}. 

The problem is complete once the velocity field satisfying the Jeans equa- 
tions under $(r) is solved for. Prom the first of these equations (cf. Binney & 
Tremaine 1987, §4.2), 



dpv, 



dr ~'~ r 



2vl 



= -p- 



d$ 

dr ' 



(16) 



for an isotropic velocity field the solution takes the form 



With a known velocity dispersion at each radius it is customary to seek closure 
by picking a form for the velocity df. We will assume a velocity dispersion 
taking locally a Maxwellian form (Hernquist 1993), 



F^{v, r) dcTv = A-n [i^j exp(-'u72v2)dv . (18) 

N-body realisations of a Maxwellian distribution would require integrating to 
infinity in velocity space. In practice bound particles all have velocities below 
the local escape velocity (= ^J—2^{r)), effectively setting an upper-limit for v 
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in (18). A proper implementation of this renormalises the velocities to preserve 
the Maxwellian character of the field, so 



where a is the one-dimensional velocity dispersion at r and — for non- 
rotating distributions. 



3.2 From spheres to axisymmetry 

To construct axisymmetric halo and bulge components, we proceeded as fol- 
lows. Our starting point is a spherical equilibrium constructed as in the above. 
With M(r) known, the potential $(r) and velocity dispersion are evaluated 
for member stars from (15) and (17). Each particle is then given a velocity 
drawn from the distribution (18). Since only one-dimensional integrals are in- 
volved, the computational requirements for this part of the calculation scales 
with oc N, where N is the total number of particles. 

We then perform a homologous transformation of the position of each particle 
to achieve the desired morphology. The case where isodensity contours are 
concentric oblate spheroids is a good starting point. The particles are then 
mapped according to {x, y) {x,y), z z\/\ — e^. The quantities <r^> and 
Tg are computed directly from M(r), prior to the mapping. 

To modify the velocity field requires evaluating the potential energy of each 
particle in the new spheroidal potential, ^^^(i?, z). The quantity V <r^> is 
computed from mass lying inside any particle radius, r, for which we de- 
fine a uniform density spheroid of axes (V <r^>, V <r'^>\J\ — e^), total mass 
M(r). Equations (8) and (9) then yield the correct potential at (i?, 2;) up to 
quadrupolar order. By virtue of the virial theorem half of the potential energy 
accrued by the stars through the coordinate transformation should be invested 
in kinetic energy. Thus 



for each particle i. The velocity vector v\ must satisfy (13), which is the case 
when 



v\r) = 3 a\T) , 



(19) 




(20) 
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In the above the primed quantities refer to the transformed, final state. The 
multiphcative factor ensures the new velocity gives the kinetic energy (20). 
The new potential could have been evaluated from (10), by computing the 
inertia tensor inside the spheroid (r, r\/l — e^), for each particle. However (9) 
and (8) require only evaluating trigonometric functions and the algorithm is 
therefore optimised for the transformation we have performed here. When 
the eccentricity of the isodensity contours vary with position, one requires 
averaged eccentricities in (13) for every mass shell, which is provided through 
(10) by the tensor of inertia. 



3.3 Numerical Setup 

The algorithm we have described was tested using the nested-grid integrator 
SUPERBOX (Fellhauer et al. 2000). The tests we conducted used in total N — 
50000 particles. Since spherical equilibria provide the work horse of the scheme, 
we first constructed a spherical equilibrium to tune up the code. 

The spherical profile used is the Dehnen family (14) with 7 = 1 and a — 1/10. 
The models were all truncated at Rcut — 1- We then evolved the distributions 
for 3 to 4 circular orbits at the half-mass radius, each revolution corresponding 
to time intervals 4tcr = 4y^37r/16(G'p)"-'^, where tcr is the dynamical time 
and p the averaged density at the half-mass radius. Setting G = 2 we have 
Atci — 0.727.. for a single revolution about the half- mass volume. At the edge 
of the system an orbit takes ~ 2.05.. units of time for a revolution. 

SUPERBOX is a time-centred leap-frog integrator. The timestep we selected 
= 0.006 X tcr] the radius where this timestep matches the local dynamical 
time for an Dehnen 7 = 1 model is at (roughly) 1 x 10~^ x a, enclosing ~ 10~^ 
percent of the mass. The dynamics is therefore well resolved in time for all 
particles in the simulations. In general the central cusp is poorly resolved by 
grid codes due to finite linear resolution. However the three levels of resolution 
allowed by SUPERBOX means that we can focus on the central region and boost 
resolution as needed. Setting the total grid size — 1/2 we then find for an 32- 
mesh grid a hnear resolution I 1/2 x 1/32 = 0.016 such that 2 percent 
of the total mass falls inside the innermost grid points. This would minimise 
structural evolution, however the coarse central density due to finite resolution 
means that velocity field and potential no longer match there, resulting in a 
small degree of relaxation at the centre. 
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The left-hand panels on Fig. 2 show the results of evolving a spherical Dehnen 
7 = 1 model. The figure displays the 10-percentile Lagrange radii of the mass 
sorted in spherical shells. We note that systematic fiuctuations are low for the 
innermost mass shells. However near the edge appreciable evolution can be 
seen, at a level of 15% for the 90% mass radius. Bearing in mind this caveat, 
we then applied the same numerical setup to flattened Dehnen models. To 
preserve the good resolution of the dynamics in time, we enforced identical 
central density for all models by an homologous transformation of the particle 
positions, so the dynamical time tcr remains the same. 



3.4 Results for axisymmetric models with 7 = 1 



Non-rotating self-gravitating equilibria of aspect ratios c/a > 1/3 are suscepti- 
ble to bending modes of instability (Combes et al. 1990, Sellwood & Wilkinson 
1993). No elliptical galaxy has a projected aspect ratio flatter than E7 (Bin- 
ney & Merrifleld 1998). Furthermore the narrow stream of debris from the 
Sagittarius dwarf suggests that the halo of the Milky Way may be as round as 
or rounder than El (Ibata & Lewis 1998). Therefore we consider a flattened 
model with c : a = 1 : 3 as a limiting case to test the algorithm. This is done 
flrst, before investigating the effects of varying aspect ratio and power index 
7- 

Figure2(b) graphs the Lagrange radii of an 1:3 flattened Dehnen model. The 
mass proflle has been binned in spherical shells to facilitate comparisons with 
the spherical case. Assessing the overall setup, we find as expected a degree of 
relaxation in the innermost mass bins. Close inspection shows fiuctuations of 
amplitude ~ 0.01 forming at all Lagrange radii. These fiuctuations propagate 
inside out, and result from stars leaving the central region. For the model at 
hand we estimate that 2.8% of the total mass gives rise to them. At the end 
of this simulation only 0.52% of stars had left the volume of the simulation, 
therefore we conclude from this exercise that fa 2% of kinetic energy is not 
assigned correctly by our application of (20) and leads to large errors in the 
binding energy of relatively few particles. 

We monitor in time the morphology of the model defining the mass- weighted 
principal axes 



where {/j} are the eigenvalues of the inertia tensor and {i.j.k} = {x.y.z} 
and permutations thereof (Goldstein 1980, §5-3; see also Kroupa 1997). Par- 
ticles are first selected in spherical Lagrangian shells, from which the /j's are 
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computed. A first estimate of tlie principal axes of the mass distribution is 
computed from (21). A second selection of particles is then made, with the 
bounding volume now assuming the ellipsoidal shape defined by the a^'s. The 
inertia tensor and eigenvalues are recomputed from this new selection to ob- 
tain the final estimate of the system axes. Figure 3 (left-hand panel) graphs 
the evolution of the quantity 203/(01 + 02) for two mass bins, the 30% and 60% 
Lagrange radii. For clarity the 30% results were divided by a factor two. We 
find the relaxed configuration giving nearly the exact same quantities after a 
half-unit of evolution. It is noteworthy that relaxation was more severe in the 
central region, where the computed aspect ratio drops from 0.18 initially to 
0.16 at later times. The right-hand panel shows the velocity field for the 60% 
innermost particles. Here there are virtually no signs of evolution or trends 
beyond the level of root-iV fiuctuations. 



3.5 Models wtth ^ 1 



We investigated the range of applicability of our approach by evolving Dehnen 
models with different power indices 7. The solution (7) for the velocity eUipsoid 
applies to harmonic potentials of uniform density. Since we opted instead 
to implement (13), we ask whether equilibria can be setup this way with a 
harmonic central core, i.e. for models with 7 = 0. 

Figure 4 (left-hand panel) illustrates the impact of the central harmonic core. 
We have graphed the same quantities as on Fig. 3, using the same numerical 
setup to validate the comparison. The oscillations of the 60% axis ratio seen 
early on in the simulation compare in magnitude to those observed for the 
7=1 case. However in the inner region the ratio varies from 0.18 initially, 
to 0.15 after one time unit of evolution, or down 17%, somewhat more than 
what was observed for the same volume in the case of 7 = 1, where variations 
were of ~ 11% relative magnitude. 

The multipolar expansion of the gravitational field suffers less from trunca- 
tion if the matter is more concentrated about the centre of gravity. Therefore 
Jaffe's (1983) profiles (Dehnen 7 = 2) provide another test since they are more 
concentrated. The right-hand panels on Fig. 4 show indeed that fluctuations 
due to dynamical relaxation are much reduced for the period t < 1/2. The 
principal axes ratio fluctuates between 0.17 ± 0.01 or 5.9% for the inner 30% 
mass radius. At larger volume it fluctuates initially by a similar relative frac- 
tion. We would expect this property to prove helpful when modelling dark 
haloes as isothermal gas, for which asymptotically p oc at large distances. 
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4 Multi-component Models 



The map (13) provides equilibria for one-component self-gravitating spheroidal 
galaxies. We now generalise to multi-component galaxies. Our goal is to verify 
that the embedding of new components does not introduce further evolution in 
the structure of the galaxy. In an effort to bring out algorithmic, as opposed 
to dynamical, biases, we have put together a bulge-disc-halo model galaxy 
where the bulge and halo have the same spherical Dehnen 7 = 1 mass profiles. 
We then opted to combine a spherical bulge, exponential disc and axisymmet- 
ric halo, requiring respectively zero-, high- and low-order corrections to the 
monopole for a precise determination of the potentials. We modified a version 
of BUILDGAL kindly provided to us by L. Hernquist and based on his 1993 
algorithm. We implemented the map of Section 2.2 for axisymmetric galactic 
potentials, taking into account the hierarchy of individual components. 

4.1 Disc, bulge & halo galaxy 

While not attempting to model the Galaxy, we set up model galaxies with 
relative masses and lengths in rough agreement with those of Freudenreich 
(1998) for the Milky Way (see Sackett 1997, Binney & Merrifield 1999, §10). 
In model units, the bulge has a mass and radius = 1; the core length Tc = 0.10. 
The exponential disc has a radial scale length h = 1.0, scale height Zo = 0.3, 
total radius = 10 and mass = 3.0. We chose scale lengths h and Zo large 
compared with the linear resolution / = 0.016 in the inner region of our grid 
code. 

The heavier spheroidal halo was given an 1:2 aspect ratio as described in 
Section 2.2. We took z as the symmetry axis so the spheroid equator lies in 
the x-y plane. The halo has a total radius = 20.00, mass = 30.00, and we 
chose a halo core length Tc = 2.0. The halo and disc masses are in the ratio 
4:1 at the edge of the disc, so the halo drives the dynamics there but not 
overwhelmingly so. This is done on purpose, with a view to perturb the halo 
equilibrium, as it must adjust to the disc gravity. 

4-2 Computing the potentials 

Experience and pragmatic considerations led to the following simplifying tricks. 
When computing the disc feedback on embedded bulge particles, we have kept 
the direct-summation algorithm from the original code, since the dimensions of 
the bulge are comparable to the disc scale length. Hence a multipole expansion 
would require high-order series expansion which we deemed not essential, as 
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the number of bulge particles is less or comparable to the number of disc par- 
ticles, which already gives t oc N]^ scaling of computational time (cf. equation 
2). We found the direct-summation computation, contributing a term oc N^Ni, 
in the total computational time, gives better result than a expanding the disc 
potential to high-order, although a crude evaluation using (8) also provides a 
sensible equilibrium - see below. The feed-back of the bulge on disc particles 
however is accounted for already when the monopole term of its potential is 
added to that of the disc. 

The halo and disc interactions require some care. We found it sufficient for an 
axisymmetric halo to expand its potential to quadrupolar term (cf. [10]) and 
adding this to the disc's (and the bulge's) self-gravitating potentials. Equation 
(10) is easily differentiated in three dimensions. The kinetic energy is then 
added to the disc and bulge particles for each degree of freedom in the same 
proportion as the components of the halo gravitational forces (to account for 
the work done in reshaping the system). In this way, relatively more kinetic 
energy is invested in the x-y plane, when compared with the z-direction, to 
fight off the enhanced gravity due to the flattened halo in the equatorial plane 
of the spheroid. As the rms disc aspect ratio is close to 1:12, the boost in 
disc kinetic energy required to maintain equilibrium was added to the circular 
motion at cylindrical radius R for each disc particle. Since the streaming 
motion at R is known, energy can be directly added to each particle once (20) 
has been computed. 

As halo particles spend httle time near or in the plane of the disc, we treated 
the disc as a highly flattened oblate spheroid, for which the potential is com- 
puted from (8). To do so we averaged over the vertical structure of the disc, 
however the original scales of height and radial decay are preserved if we use 
as spheroid axes the mass- weighted disc parameters <R> and <z>. This 
effectively removes the term oc N^Nh of the total computational time (2). 



4-3 Results 

For comparison purposes we display two models on Fig. 5. The first one has 
both spherical bulge and halo, while the second has an oblate spheroidal halo 
of aspect ratio 1:2. The motivation for doing this is clear: for spherically sym- 
metric components, the multipole approach becomes exact, since only the 
monopole term contributes to the potential of each component. Models of a 
disc galaxy with spherical bulge and halo thus minimise relaxation effects. 

The structure of the bulge is well accounted for by sorting the mass on spher- 
ical Lagrangian shells. For the case of the disc, of mean aspect ratio ~ 1 : 12, 
the Lagrangian shells effectively follow the cyhndrical radii. Contrasting the 
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two models, with and without spherically symmetric halo, we find little ef- 
fect on the radial structure of disc and bulge introduced by our approximate 
treatment. In both the cases, radial oscillations are seen to propagate from the 
inside, outwards, while taking several dynamical times to dampen. On the fig- 
ure the arrows indicate the dynamical time ^cr at the half-mass radius (bulge 
and halo) or at R — h (disc) of individual components in isolation. These 
radial oscillations give rise to transient patterns on visual inspection. We note 
that after a few revolutions of the disc at R = h, much of the fluctuations 
have disappeared or are much attenuated. 



The moments of inertia give insight into the structure of individual galaxy 
components. On Fig. 6 we display the moments of inertia, principal axes and 
velocity dispersion components of the bulge and halo for the case where the 
halo is an oblate spheroid. We find as on Fig. 2 slight changes in the ratio of 
principal axes of the halo. However, after settling, the structure remains close 
to the one specified on input. The bulge, on the other hand, is clearly stretched 
out of spherical symmetry, due to the combined pull of the disc and halo: at 
the edge of the bulge, R = 1, the halo mass ~ 10/3 is already several times 
that of the bulge, which is not fully self-gravitating anymore. The aspect ratio 
of the bulge in equilibrium is close to 0.65/0.75 = 0.87, when we would want 
unity. The external forces acting on the bulge leads to expansion, as seen from 
the increased moments of inertia (top left-hand side panel). A quick solution 
to counter this efi'ect would be to stretch the bulge initially to mildly prolate 
structure of axis ratio a : c ~ 1 : 1.1, so as to relax to near sphericity. We have 
not resorted to this therapy and chose instead to illustrate the limitations of 
the algorithm. 



Turning to the vertical structure of the disc, we averaged the height \z\ of 
stars inside R — 2h, since this region is well resolved by our grid code and 
errors due to numerical integration are reduced. We computed both average 
vertical height and the one-dimensional velocity dispersion cTz for ten snapshots 
covering ten time units of evolution, or four revolutions. Table 1 gives the 
results. Early in the evolution of the disc, the vertical structure shrunk by 
15%, as the mean height drops from 0.2 to ^ 0.17. We observed a decrease 
of the standard deviation as well, of the same magnitude, which leads us to 
conclude in a global contraction of the vertical disc structure. This is matched 
by a rise in velocity dispersion of the same relative magnitude (see Table 1). 
Since < >~ and is much less than the vertical velocity dispersion, the 
bulk motion remains negligible thoughout. 
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4-4 Benchmarks 



Approximate methods such as the one we have presented would not be of 
interest without a handsome pay-out at the computer. We find the scahng 
properties of the method with a series of numerical tests. Specifically, we first 
bench-marked the original version of the code BUILDGALfor spherical bulge 
-|- disc galaxy models. We then constructed an axisymmetric bulge using the 
full algorithm, which computes the response of the disc to bulge potential 
exactly. Since bulge and halo are constructed in the same way in buildgal , 
and the disc embedding done in the same exact fashion, we may relabel the 
bulge as halo. This is so when referring to BUILDGAL below. In a second set of 
tests, we repeated these experiments, using the new version of the code. The 
same default disc parameters were used in all our tests, and haloes were given 
Hernquist (1990) mass profiles with Tc = 0.1, and truncated at the disc edge. 

All benchmarks refer to an 333 MHz G3 Apple processor operating under 
S.u.S.E. Linux 6.4 OS; we compiled buildgal with the widely available GNU 
Fortran compiler g77 -04. The performance were done with disc and bulge 
particle numbers, N^^Ni,, in the range 5,000 to 5 x 10^ (disc); and 5000 to 
A^b = 5 X 10^ (bulge or halo). 

We solved for the operation times t^^i in (2), first as a set of linear equations 
by constructing stand-alone disc and bulge (no feedback terms present), then 
re-doing the calculations with disc and bulge together. Our results are hsted 
in Tables 2 and 3. The most important point to notice is the large increase in 
computer time for models with an axisymmetric bulge, which scales effectively 
with the square of particle number. The turn-over to quadratic dependence is 
approximately 2 10'^ particles. The validity of (2) is inferred from the conver- 
gence of the operation times t^/s with increasing N. We have extrapolated to 
infinite particle numbers using the polynomial fit routine pzextr of Press et 
al. (1992). The numbers are truncated to the last significant digit according to 
error estimates obtained from the results of the two highest-order polynomial 
fits. 

Repeating this exercise with the new algorithm, we concentrated on the case 
of an axisymmetric halo -|- disc, since the case of a spherical halo would give 
essentially the same tx/s. We therefore have to compute only the extra com- 
putational time required for the transformation of the halo to axisymmetry 
and its coupling with the disc. We find as expected that the time required 
to adjust the disc and halo velocity fields scale in proportion to their respec- 
tive particle numbers. A comparison of total computational time shows that 
axisymmetric compound models require only a modest increase in computer 
time when compared with spherically symmetric components. The new code, 
MaGalie (for Mache Galaxie), can thus construct galaxies with bulge, disc 
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and halo, all in axisymmetry, at virtually no extra costs. 



5 CONCLUSIONS 

We have presented a method of constructing approximate equilibria in good 
agreement with those sought. In particular the algorithm scales linearly with 
bulge and halo particle numbers, when these are axisymmetric or spherical. 
The approach takes into account implicitly the relative scales of components 
involved, and treats the interactions between components by expanding the 
potential to required order (Fig. 1). 

The errors introduced from treating non-spherically symmetric components 
was found to be no worse than the case where only spherical components 
surround a disc. The required computational time for this algorithm scales 
with particle numbers as 



^cpu — tden,h + -^/i ivel,h + -^fe icar,h + -^d ^den.d + -^d 'tvel,d + tcor,d (22) 

where the parameters tx,i are given in Table 3 for a G3 processor (see Sec- 
tion 4.4). With this algorithm non-spherical haloes surrounding galaxies can 
be constructed as a matter of routine. The results of a parameter survey of 
orbital decay in aspherical haloes will be presented in a separate contribution 
(Peiiarrubia et al. 2000). 

We have limited our discussion to prolate spheroidal potentials. The approach 
is readily applicable to prolate or triaxial structures however. To construct 
triaxial equilibria, one would map the axes of the velocity ellipsoid as in (13) 
to the geometry of the triaxial body. Equations (8) and (10) are easily gener- 
alised to triaxial systems (BT-l-87). Likewise, net rotation of individual com- 
ponents is implemented by alignment of particle angular momenta (see e.g. 
Lynden-Bell 1962, Hernquist 1993). Models of barred galaxies may possibly be 
constructed by converting random into bulk motion of a triaxial component, 
but we have not implemented this. 

The disc self-gravity now represents the bottleneck through the N]^ term in 
equation (22). For situations where the disc substructure can be neglected, 
this dependence may be removed by constructing a disc from flattening a 
spheroid to small aspect ratio, then converting random motion to provide ro- 
tational support. When the disc structure matters to the dynamics however, 
this approach is unlikely to yield sensible disc equilibria, because the vertical 
and radial motions in thin structures require fine-tuning the velocities locally, 
hardly in line with such a brute-force approach. Note however that one may 
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construct a disc similar to the Milky Way's by considering the thin and thick 
discs as separate components. These may then be combined in the fashion 
described in Section 2. Sellwood & Wilkinson (1993) discuss methods of con- 
structing stable thin galactic discs comprehensively. 
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Note added in proof - After this article was submitted, we found a natural, 
quick-fix solution to the Wj' scaling of the thin disc construction algorithm 
(cf. Tables 2 and 3): it consists in adding together m realisations of an n- 
particle disc to reach the desired total particle number Na = m x n. The total 
cpu then scales with m x ri^ n x Nd, or linearly with Nd- Though n and 
hence the coefficient of proportionality should be reasonably large to dampen 
numerical noise, the approach is well-suited to parallel algorithms. 
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Table 1 

Time-evolution of the vertical structure of an embedded disc in a flattened spheroid. 
Only particles inside cylindrical radius R = 2h, where h is the radial disc scale 
length, are included when computing mean height and velocity dispersion. The 
standard deviation 5z = <(|-2|— <\z\>Y>. 
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Component 3 



Component! 



Fig. 1. Coupling galaxy components of dif- 
ferent linear sizes and masses. A low ac- 
curacy treatment of the total potential 
is sufficient to maintain the largest, most 
massive component in equilibrium, while 
it induces a relatively stronger response 
from the smaller components, requiring 
higher-order treatment of the potential for 
these components. The strength of the 
feedback is shown here as light and heavy 
arrows. 
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(A) Total CpU — Nh X {tden,h + tvel,h) + 2 Nd Nfi tcor,h + -^d tden,d + -^J ^t;e/,(i 



Nh 


tden,h 


tyeljh 


tcor,h 


Nd 


tden,d 


tvel,d 


CPU [seconds] 




[10-5 s] 


[10-5 s] 


[10-5 s] 




[10-5 s] 


[10-5 s] 


True 


Computed 


510^ 


3.80 


6.70 


0.096 




27.2 


0.368 


141.9 


83.7 


10^ 


3.90 


6.10 


0.092 


510^ 


27.2 


0.368 


186.4 


141.1 


lO'' 


1 


1 


0.096 


10" 


27.3 


0.222 


417.7 


231.4 


10^ 






0.106 


210-^ 


27.4 


0.159 


1066.5 


862.0 


10 


3.90 


6.10 


0.107 


5 10 


27.3 


0.112 


3884.6 


3653.9 


210^ 


3.80 


4.90 




105 


26.8 


0.102 


10228.5 


10027.9 


510^ 


3.80 


4.10 










3.95 


3.64 


10^ 


3.80 


3.80 










7.60 


7.28 


510^ 


3.70 


3.60 










36.50 


36.40 


oo 


3.64 


3.54 


0.114 


oo 


26.5 


0.100 






(B) Total cpu = 


Nh tden,h 


+ N^ tvel,h + '^NdNh tcor,h 


+ Nd tden,d 


+ Nj tyel,d 




Nh 


tden,h 


tyeljh 


tcor,h 


Nd 


tden,d 


tvel,d 


CPU [seconds] 




[10-5 s] 


[10-5 s] 


[10-5 s] 




[10-5 s] 


[10-5 s] 


True 


Computed 


510=^ 


12.4 


(0.120) 


0.096 


510^ 


27.2 


0.368 


172.0 


110.4 


10^ 


12.4 


0.115 


0.092 


510^ 


27.2 


0.368 


301.6 


247.6 


10^ 




1 


0.096 


10" 


27.3 


0.222 


532.9 


437.9 


10^ 






0.106 


210" 


27.4 


0.159 


1181.7 


968.5 


10^ 


12.4 


0.115 


0.107 


510" 


27.3 


0.112 


3999.9 


3760.5 


210'' 


12.5 


0.108 




105 


26.8 


0.102 


10661.3 


10453.0 


210^ 


12.5 


0.108 










434.5 


426.5 


510^ 


12.4 


0.107 










2681.2 


2656.9 


oo 


12.4 


0.106 


0.114 


oo 


26.5 


0.100 







Tabic 2 

Benchmarks of the code BUILDGAL for (a) spherical halo + disc and (b) axisymmetric 
halo + disc models. Notation refers to (2); tcor,d = tcor,h by symmetry and has been 
omitted in the table. We computed CPU values from (2) using the asymptotic limits 
shown in bold face; true CPU corresponds to clocked time. 



23 



Total CpU = Nh tden,h + Nh tyel,h + Nh tcor,h + -^d tden,d + tyel,d + tcor,d 
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Table 3 

Benchmarks of the new code MaGalie for spherical or axisymmetric halo + disc 
models. Notation refers to (2). 
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a) Spherical 



b) c:a - 1 :3 




0.0 1.0 2.0 3.0 4.0 0.0 1.0 2.0 3.0 4.0 

time [m.u.] time [m.u.] 



Fig. 2. (a) Evolution of a spherical Dehnen 7 = 1 model, (b) Same model, but 
flattened to achieve an aspect ratio = 1:3. The 10-percentile Lagrange radii are 
displayed versus time, where one unit of time corresponds to a full revolution at 
half-mass. 
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Fig. 3. Left panel: structural evolution of a flattened 1:3 Dehnen 7=1 model. The 

ratio of minor to major axes are computed at two different Lagrange radii enclosing 
30% or 60% of the total mass. Note: the results at 30% Lagrange radius were divided 
by two to avoid overlap with the others. Right panel: components of the velocity 
dispersion tensor (of. Eq. 21) for the innermost 60% particles. The lower curve is 
the z-component. The velocity dispersions were computed along the eigenvectors of 
the rotational inertia tensor. 
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Fig. 4. As for Fig.3. Time evolution of two flattened Dehnen models of aspect ratio 
1:3, with 7 = (left) and 7 = 2 (right). The principal axes (21) were computed at 
the 30% and 60% spherical Lagrange radii. The 30%-results were divided by two to 
distinguish the two curves on each panels. 
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a) Spherical Dehnen y = 1.0 bulge ; r^, = 1/10, N,, = 10,000 



a) Spherical Dehnen y = 1.0 bulge ; r^, = 1/10, N,, = 10,000 
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b) Exponential disc ; h = 1, = 10,000 
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b) Exponential disc ; h = 1, = 10,000 
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c) Spherical halo ; = 2.0, = 50,000 



c) Spheroidal 1:2 halo ; r^ = 2.0, = 50,000 
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Fig. 5. Left panels: evolution of the Lagrange radii (including 10, 20, 30 ... percent 
of the component's mass) for a three-component model with spherical bulge and 
halo. Right panels: as on the left, but with a spheroidal halo of aspect ratio 1:2. The 
arrows indicate the circular period at the 50% Lagrange radius (bulge and halo) 
and at the radial scale length h (disc, middle panel). 
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Spherical Dehnen y = 1.0 bulge ; r^, = 1/10, = 10 
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Spheroidal 1 :2 halo; r^, = 2, Nf, = 5x1 
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Fig. 6. Time-evolution of the morphology of individual components. Left-hand pan- 
els: central bulge parameters; right-hand panels: halo parameters. The quantities 
were measured at the 50% spherical Lagrange radius. The principal axes were com- 
puted as in (Fig. 3). The XYZ Inertia are the eigenvalues of the rotational inertia 
tensor (Goldstein 1980). 
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